%%%%%%%%%%%%%%%%%%
% Summary stats and results from running return factor model
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%%

minN=200;winth=0.01; nx=5; nlag=12;


dds=Z.dds;
ids=Z.ids;
dd=Z.dd; T=length(dd);
IV=mean(Z.IVV(:,1:2,1),2); %ATMV
y=mean(Z.NRV(:,1:2,1,3),2)*100; %straddle return, daily updated hedged
x=Z.Factors;%[hedging cost, vol risk, jump];
HVV=Z.HVV; %historical vol at 3 horizons


Bfv=NaN(T,nx);Afv=NaN(T,1);R2v=NaN(T,1);
statv=NaN(11,nx+2,T);
sx=NaN(T,nx);bi=NaN(T,4);
nmd=5;qx=NaN(nx,nmd,T); qy=qx;qv=qy;VB=NaN(T,nx);
parfor t=nlag+1:T
    [Bi,nobt,Af,Bf,R2s,stats,sxt,sxt2,bit,r2it,qxt,qyt,qvt,vbt] =FunDailyORSestimation(t,y,x,HVV,IV,ids,dds,dd,nlag,nx,winth,minN);
   
    Bfv(t,:)=Bf; R2v(t)=R2s;
    Afv(t)=Af;
    statv(:,:,t)=stats;
    sx(t,:)=sxt;
    bi(t,:)=[bit;r2it];
    qx(:,:,t)=qxt; qy(:,:,t)=qyt;qv(:,:,t)=qvt; VB(t,:)=vbt;
end

fl=[Afv,Bfv,R2v*100];

Tf=sum(isfinite(fl));
flstats=[nanmean(fl)', nanstd(fl,0)',sqrt(Tf)'.*nanmean(fl)'./nanneweystd(fl)',sqrt(12).*nanmean(fl)'./nanstd(fl)', nanacf(fl,1)', skewness(fl,0)', kurtosis(fl,0)'-3] ;

%intertemporal relation
ax=reshape(statv(1,2:4,:),3,T)';
ii=isfinite(sum([ax,fl],2)); Ts=sum(ii);
yx=[fl(ii,2:4)];%returns on 3 factor portfolios
xx=ax(ii,:); %average levels of 3 risk factors only
xx=sx(ii,1:3); %standard deviaiton of the 3 risk factors.
xs=[ones(Ts,1),(xx-repmat(nanmean(xx),Ts,1))./repmat(nanstd(xx),Ts,1)];
Bs=xs\yx;
es=yx-xs*Bs;
sbs=sqrt(diag(inv(xs'*xs))).*neweystd(es);
tbs=Bs./sbs;
atb=abs(tbs);

AR2=(1- sum(es.^2)./(Ts-4)./std(yx).^2);
a=[];
for j=2:4
    a=[a,Bs(j,:)',atb(j,:)'];
end
aav=a;

disp('Table 2, Summary Stats of risk exposures on delta-hedged option investments');
Astats=nanmean(statv,3)';

xname={'Delta hedging cost','Stochastic volatility','Random jumps','Historical premium','Volatility premium'}';
fprintf(1, '  \\\\ \n');
for j=1:nx;fprintf(1, '%40s ', xname{j});fprintf(1, [repmat(' & %8.2f  ',1,2),'&',repmat(' & %8.2f  ',1,5),'&',repmat(' & %8.2f  ',1,2) '  \\\\ \n'], Astats(j+1,(1:end-2)));end
fprintf(1, '   \n');

disp('Table 3, Delta hedged option excess returns and return volatilities across risk quintiles');
mqy=nanmean(qy,3);
hl=reshape(qy(:,end,:)-qy(:,1,:),nx,T)';
mhl=nanmean(hl);
shl=nanstd(hl); Ti=sum(isfinite(hl));
thl=sqrt(Ti).*mhl./shl;

disp('A. Option mean excess returns');
fprintf(1, '  \\\\ \n');
for j=1:nx;fprintf(1, '%40s ', xname{j});fprintf(1, [repmat(' & %8.2f  ',1,5),' && %8.2f & ( %8.2f )  \\\\ \n'], [mqy(j,:),mhl(j), thl(j)]);end

disp('B. Option return volatilities');
mqv=nanmean(qv,3);
hlv=reshape(qv(:,end,:)-qv(:,1,:),nx,T)';
mhlv=nanmean(hlv);
shlv=nanstd(hlv); Ti=sum(isfinite(hlv));
thlv=sqrt(Ti).*mhlv./shlv;

fprintf(1, '  \\\\ \n');
for j=1:nx;fprintf(1, '%40s ', xname{j});fprintf(1, [repmat(' & %8.2f  ',1,5),' && %8.2f & ( %8.2f )  \\\\ \n'], [mqv(j,:),mhlv(j), thlv(j)]);end
fprintf(1, '   \n');


disp('Table 4. The risk-return relation on delta-hedged option investments');
ax=reshape(statv(1,1,:),T,1); Tf=sum(isfinite(ax));
axs=[nanmean(ax), nanstd(ax,0)',sqrt(Tf)'.*nanmean(ax)'./nanneweystd(ax)',sqrt(12).*nanmean(ax)'./nanstd(ax)', nanacf(ax,1)', skewness(ax,0)', kurtosis(ax,0)'-3] ;
fprintf(1, '  \\\\ \n');
fprintf(1, '%40s ', 'Average');fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,5),' \\\\ \n'], axs);

xnamel={'Intercept','Hedging cost','Volatility risk','Jump risk','Historical premium','Volatility premium','Adjusted $R^2$'}';
mVB=[NaN;100*nanmean(VB)'; NaN];
fprintf(1, '  \\\\ \n');
for j=1:1;fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,5), ' \\\\ \n'], [flstats(j,:)]);end
for j=2:nx+1;fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,6), ' \\\\ \n'], [flstats(j,:),mVB(j)]);end
for j=nx+2:nx+2;fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,5), ' \\\\ \n'], [flstats(j,:)]);end
fprintf(1, '   \n');

disp('%Table 5: intertemporal relation');
pname={'Hedging cost','Volatility risk','Jump risk'}';nps=length(pname);
fprintf(1,' \n')
for j=1:3
    fprintf(1, '%30s ', pname{j});fprintf(1, [' & ', repmat('  %8.2f  & ( %8.2f ) && ',1,nx-3),'  %8.2f  & ( %8.2f )  && %8.2f  \\\\ \n'], [aav(j,:),100*AR2(j)]);
end



